Permutation Function

The Function

We create an R version of Algorithm 7.6.1 from Thisted’s draft. This version attempts to use a for loop instead of a repeat loop.

  x <- 1:10

  ### Randomly Permutate a vector
  RanShuffle <- function(x){
   ###initiate random variable
  
    for(i in 1:length(x)){
    ###Pick A Random Index From 1 to i 
      randomInd <- floor(runif(1,1,i))
    ###Swap X[i] Elelment with element at random index
      temp<-x[i]
      x[i] <- x[randomInd]
      x[randomInd]<-temp
    } 
    return(x)
    
  }

  RanShuffle(x)

##  [1]  2  6  9  8  1  7 10  5  4  3

  RanShuffle(x)

##  [1]  4  1  8  9  3  7  5 10  6  2

  RanShuffle(x)

##  [1]  5 10  8  9  6  4  3  1  2  7

Check the Function

We look to see if it works.

  ### Create a function that evaluates the z values as factors so that
  ### all values listed in the input variable, y, are used
  table.factor <- function(z, lvls=y){
                                       table(factor(z, levels=lvls))
                                     }

  n <- 10
  nreps <- 10000
  y <- 1:n
  z <- RanShuffle(y)
  for (i in 2:nreps){
    z <- rbind(z, RanShuffle(y))
  }
  z_counts <- apply(z, 2, table.factor)
  z_counts

##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## 1     0 1149 1142 1116 1059 1023 1106 1153 1164  1088
## 2  1106    0 1120 1108 1132 1090 1112 1123 1093  1116
## 3  1066 1119    0 1126 1074 1145 1115 1117 1128  1110
## 4  1084 1131 1140    0 1101 1087 1099 1144 1054  1160
## 5  1113 1185 1112 1140    0 1093 1091 1127 1045  1094
## 6  1152 1058 1120 1108 1158    0 1073 1069 1125  1137
## 7  1136 1080 1099 1083 1116 1162    0 1077 1169  1078
## 8  1133 1122 1086 1108 1138 1109 1086    0 1142  1076
## 9  1116 1097 1091 1083 1122 1140 1157 1053    0  1141
## 10 1094 1059 1090 1128 1100 1151 1161 1137 1080     0

A “Heatmap”

Maybe a heatmap will help us see if there are problems.

  p_load(reshape2, ggplot2)
  melted_z_counts <- melt(z_counts)
  head(melted_z_counts)

##   Var1 Var2 value
## 1    1    1     0
## 2    2    1  1106
## 3    3    1  1066
## 4    4    1  1084
## 5    5    1  1113
## 6    6    1  1152

  names(melted_z_counts) <- c("Location","Value","Count")
  
  ### Plot the deviations of the observed (count) away from the 
  ### predicted (nreps/n)
  ggplot(data = melted_z_counts, aes(x=Location, y=Value, fill=Count-nreps/n)) + 
  geom_tile() + scale_fill_gradient(low = "red", high = "green")

Clearly, we have a problem.